This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
library(Seurat)
## Loading required package: SeuratObject
## Loading required package: sp
##
## Attaching package: 'SeuratObject'
## The following objects are masked from 'package:base':
##
## intersect, t
library(ggplot2)
library(RColorBrewer)
load("/Users/w435u/Documents/ST_SC/DATA_STSC_CAO/Seurat_C2.RData")
SpatialFeaturePlot(st_dataset, "pred_st")+
scale_fill_gradientn(colours = rev(brewer.pal(n = 11, name = "RdBu"))[1:11])
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
You can also embed plots, for example:
df <- st_dataset@meta.data
# Create plot
library(ggpubr)
p1 <- ggplot(df, aes(x = cnv_score, y = pred_st)) +
geom_point(alpha=0.5) + # Scatter plot
geom_smooth(method = "lm", se = TRUE, color = "#ae0000") + # Fitted regression line
stat_cor(method = "pearson", digits = 2) + # Correlation and p-value
theme_minimal()+
scale_y_continuous(limits = c(0, 1))+ylab("Malignancy") + xlab("CNV score")
p2 <- ggplot(df, aes(x = tumor_sig1, y = pred_st)) +
geom_point(alpha=0.5) + # Scatter plot
geom_smooth(method = "lm", se = TRUE, color = "#ae0000") + # Fitted regression line
stat_cor(method = "pearson") + # Correlation and p-value
theme_minimal()+
scale_y_continuous(limits = c(0, 1))+ylab("Malignancy")+ xlab("Tumor signature")
p3 <- ggplot(df, aes(x = pEMT1, y = pred_st)) +
geom_point(alpha=0.5) + # Scatter plot
geom_smooth(method = "lm", se = TRUE, color = "#ae0000") + # Fitted regression line
stat_cor(method = "pearson") + # Correlation and p-value
theme_minimal()+
scale_y_continuous(limits = c(0, 1))+ylab("Malignancy")+ xlab("pEMT score")
p <- p1+p2+p3
pdf(file=paste0("~/Documents/GitHub/SCTP_iMETA/Fig2/corr_", "Q3", ".pdf"), width=12, height=6)
print(p)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 364 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 364 rows containing non-finite outside the scale range
## (`stat_cor()`).
## Warning: Removed 364 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
dev.off()
## quartz_off_screen
## 2
p
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 364 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 364 rows containing non-finite outside the scale range
## (`stat_cor()`).
## Warning: Removed 364 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Dot
plot of different
pred_st <- st_dataset@meta.data$pred_st
pred_cat <- as.character(cut(pred_st,
breaks = c(-Inf, quantile(pred_st, c(0.25, 0.5, 0.75, 1))),
labels = c("C1", "C2", "C3", "C4")))
st_dataset <- AddMetaData(st_dataset, metadata = pred_cat, col.name = "predcat")
#tumor_enriched <- FindMarkers(st_dataset, ident.1 = "C4",ident.2 = "C1", group.by = 'predcat', verbose = FALSE)
#write.csv(tumor_enriched, file="~/Documents/GitHub/SCTP_iMETA/Fig2/DEGs_P3.csv")
tumor_enriched <- read.csv( file="~/Documents/GitHub/SCTP_iMETA/Fig2/DEGs_P3.csv", row.names = 1)
genes <- rownames(tumor_enriched)[1:40]
p <- DotPlot(object = st_dataset, features =genes, group.by = "predcat")
## Warning: Scaling data with a low number of groups may produce misleading
## results
df <- p$data
#df$features.plot <- rownames(df)
exp_mat<-df %>%
dplyr::select(-pct.exp, -avg.exp) %>%
tidyr::pivot_wider(names_from = id, values_from = avg.exp.scaled) %>%
as.data.frame()
exp_mat <- exp_mat[!is.na(exp_mat$features.plot), ]
row.names(exp_mat) <- exp_mat$features.plot
exp_mat <- exp_mat[,-1] %>% data.matrix()
percent_mat<-df %>%
dplyr::select(-avg.exp, -avg.exp.scaled) %>%
tidyr::pivot_wider(names_from = id, values_from = pct.exp) %>%
as.data.frame()
percent_mat <- percent_mat[!is.na(percent_mat$features.plot), ]
row.names(percent_mat) <- percent_mat$features.plot
percent_mat <- percent_mat[,-1] %>% data.matrix()
library(viridis)
## Loading required package: viridisLite
col_fun = circlize::colorRamp2(c(-1, 0, 2), plasma(20)[c(1,10, 20)])
cell_fun = function(j, i, x, y, w, h, fill){
grid.rect(x = x, y = y, width = w, height = h,
gp = gpar(col = NA, fill = NA))
grid.circle(x=x,y=y,r= percent_mat[i, j]/100 * min(unit.c(w, h)),
gp = gpar(fill = col_fun(exp_mat[i, j]), col = NA))}
write.csv(exp_mat, file="~/Documents/GitHub/SCTP_iMETA/Fig2/Expression_FIG4H.csv")
library(ComplexHeatmap)
## Loading required package: grid
## ========================================
## ComplexHeatmap version 2.20.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
##
## If you use it in published research, please cite either one:
## - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
## - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
## genomic data. Bioinformatics 2016.
##
##
## The new InteractiveComplexHeatmap package can directly export static
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
dot <- Heatmap(exp_mat,
heatmap_legend_param=list(title="expression"),
column_title = "clustered dotplot",
col=col_fun,
rect_gp = gpar(type = "none"),
cell_fun = cell_fun,
row_names_gp = gpar(fontsize = 10),
row_km = 2,
border = "black",
cluster_columns = FALSE,column_order = c("C4","C3","C2", "C1"))
dot
pdf("~/Documents/GitHub/SCTP_iMETA/Fig2/Dotplot_Q3_degs.pdf")
print(dot)
dev.off()
## quartz_off_screen
## 2
Differential expression between different annotated groups
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
load("/Users/w435u/Documents/ST_SC/OUTPUT_CAO_new/Seurat_C4_subtumor.RData")
SpatialDimPlot(st_dataset, "T12")
#tumor_enriched1 <- FindMarkers(st_dataset, ident.1 = "T2",ident.2 = "T1", group.by = 'T12', verbose = FALSE)
#tumor_enriched2 <- FindMarkers(st_dataset, ident.1 = "T1",ident.2 = "Other", group.by = 'T12', verbose = FALSE)
#tumor_enriched3 <- FindMarkers(st_dataset, ident.1 = "T2",ident.2 = "Other", group.by = 'T12', verbose = FALSE)
#write.csv(tumor_enriched1, file="~/Documents/GitHub/SCTP_iMETA/Fig2/DEGs_Q4_1.csv")
#write.csv(tumor_enriched2, file="~/Documents/GitHub/SCTP_iMETA/Fig2/DEGs_Q4_2.csv")
#write.csv(tumor_enriched3, file="~/Documents/GitHub/SCTP_iMETA/Fig2/DEGs_Q4_3.csv")
tumor_enriched1 <- read.csv("~/Documents/GitHub/SCTP_iMETA/Fig2/DEGs_Q4_1.csv", row.names = 1)
tumor_enriched2 <- read.csv("~/Documents/GitHub/SCTP_iMETA/Fig2/DEGs_Q4_2.csv", row.names = 1)
tumor_enriched3 <- read.csv("~/Documents/GitHub/SCTP_iMETA/Fig2/DEGs_Q4_3.csv", row.names = 1)
subgenes1 <- tumor_enriched1%>% filter(p_val<10e-14)%>% filter(abs(avg_log2FC)>1.5)
subgenes1$Group <- "T2 v.s T1"
subgenes2 <- tumor_enriched2%>% filter(p_val<10e-14)%>% filter(abs(avg_log2FC)>1.5)
subgenes2$Group <- "T1 v.s Other"
subgenes3 <- tumor_enriched3%>% filter(p_val<10e-14)%>% filter(abs(avg_log2FC)>1.5)
subgenes3$Group <- "T2 v.s Other"
suball <- rbind(subgenes1, subgenes2, subgenes3)
suball$gene <- rownames(suball)
suball$label <- ifelse(suball$avg_log2FC < (-0.5), '#0068B2',
ifelse(suball$avg_log2FC > 0.5, '#B2182B',
'black'))
dbar <- suball %>% group_by(Group) %>% summarise_all(list(min=min, max=max))%>%
dplyr::select(Group, avg_log2FC_min, avg_log2FC_max)
my_col3 = c( "black","#e60012", "#00a0e9")
library(ggrepel)
g <- ggplot()+
geom_col(data = dbar, # 绘制负向背景柱状图
mapping = aes(x = Group,y = avg_log2FC_min),
fill = "#dcdcdc",alpha = 0.6, width = 0.7) +
geom_col(data = dbar, # 绘制正向背景柱状图
mapping = aes(x = Group,y = avg_log2FC_max),
fill = "#dcdcdc",alpha = 0.6, width = 0.7) +
geom_jitter(data = suball, # 绘制所有数据点
aes(x = Group, y = avg_log2FC, color = label),
size = 0.85,
width =0.3) +
geom_tile(data = suball, # 绘制中心分组标记图
aes(x = Group,
y = 0,
fill = Group),
height=0.5,
color = "black",
alpha = 0.4,
show.legend = F) +
ggsci::scale_fill_npg(alpha=0.5) + # 自定义颜色
scale_color_manual(values=rev(my_col3))+# 自定义颜色
geom_text_repel(data = suball, # 这里的filter很关键,筛选你想要标记的基因
aes(x = Group, y = avg_log2FC, label = gene),
size = 2,
max.overlaps = getOption("ggrepel.max.overlaps", default = 15),
color = 'black',
force = 1.2,
arrow = arrow(length = unit(0.008, "npc"),
type = "open", ends = "last"))+
labs(x="Cluster", y="Average logFC") +
geom_text(data=suball, # 绘制中心分组标记图文本注释
aes(x=Group,
y=0,
label=Group),
size = 3,
color ="white") +
theme_minimal() +
theme(axis.title = element_text(size = 13,color = "black",face = "bold"),
axis.line.y = element_line(color = "black",size = 1.2),
axis.line.x = element_blank(),
axis.text.x = element_blank(),
panel.grid = element_blank(),
legend.position = "top",
legend.direction = "vertical",
legend.justification = c(1,0),
legend.text = element_text(size = 13))
## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
g
## Warning: ggrepel: 783 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
pdf("~/Documents/GitHub/SCTP_iMETA/Fig2/Dotplot_Q3_degs.pdf")
print(g)
## Warning: ggrepel: 775 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
dev.off()
## quartz_off_screen
## 2
Gene expression of top genes
load("/Users/w435u/Documents/ST_SC/OUTPUT_CAO_new/Seurat_C4_subtumor.RData")
tumor_enriched1 <- read.csv("~/Documents/GitHub/SCTP_iMETA/Fig2/DEGs_Q4_1.csv", row.names = 1)
tumor_enriched2 <- read.csv("~/Documents/GitHub/SCTP_iMETA/Fig2/DEGs_Q4_2.csv", row.names = 1)
tumor_enriched3 <- read.csv("~/Documents/GitHub/SCTP_iMETA/Fig2/DEGs_Q4_3.csv", row.names = 1)
genes <- unique(c(rownames(tumor_enriched1)[1:25],rownames(tumor_enriched2)[1:25]))
p <- DotPlot(object = st_dataset, features =genes, group.by = "T12")
## Warning: Scaling data with a low number of groups may produce misleading
## results
df <- p$data
exp_mat<-df %>%
dplyr::select(-pct.exp, -avg.exp) %>%
tidyr::pivot_wider(names_from = id, values_from = avg.exp.scaled) %>%
as.data.frame()
exp_mat <- exp_mat[!is.na(exp_mat$features.plot), ]
row.names(exp_mat) <- exp_mat$features.plot
exp_mat <- exp_mat[,-1] %>% data.matrix()
percent_mat<-df %>%
dplyr::select(-avg.exp, -avg.exp.scaled) %>%
tidyr::pivot_wider(names_from = id, values_from = pct.exp) %>%
as.data.frame()
percent_mat <- percent_mat[!is.na(percent_mat$features.plot), ]
row.names(percent_mat) <- percent_mat$features.plot
percent_mat <- percent_mat[,-1] %>% data.matrix()
library(viridis)
col_fun = circlize::colorRamp2(c(-1, 0, 2), plasma(20)[c(1,10, 20)])
cell_fun = function(j, i, x, y, w, h, fill){
grid.rect(x = x, y = y, width = w, height = h,
gp = gpar(col = NA, fill = NA))
grid.circle(x=x,y=y,r= percent_mat[i, j]/100 * min(unit.c(w, h)),
gp = gpar(fill = col_fun(exp_mat[i, j]), col = NA))}
#
# write.csv(exp_mat, file="~/Documents/GitHub/SCTP_iMETA/Fig2/Expression_FIG5F.csv")
exp_mat <- read.csv("~/Documents/GitHub/SCTP_iMETA/Fig2/Expression_FIG5F.csv", row.names = 1)
dot <- Heatmap(exp_mat,
heatmap_legend_param=list(title="expression"),
column_title = "clustered dotplot",
col=col_fun,
rect_gp = gpar(type = "none"),
cell_fun = cell_fun,
row_names_gp = gpar(fontsize = 10),
row_km = 3,
border = "black",
cluster_columns = FALSE,column_order = c("T2","T1","Other"))
## Warning: The input is a data frame-like object, convert it to a matrix.
dot
pdf("~/Documents/GitHub/SCTP_iMETA/Fig2/Dotplot_Q4_degs.pdf")
dot
dev.off()
## quartz_off_screen
## 2